CopyKat——基于高通量单细胞测序方法鉴定肿瘤细胞拷贝数变异和亚克隆结构

关注生信阿拉丁，学习不打烊～

今天来看CopyKat

基于高通量单细胞测序方法

鉴定肿瘤细胞拷贝数变异和亚克隆结构

背  景

单细胞转录组技术在多种肿瘤研究中的应用，使其日趋成为解析肿瘤微环境中肿瘤细胞亚群与正常细胞亚群基因表达模式强有力的技术。

多种高通量测序技术（microdroplet systems：Drop-Seq、Indrop5、10X Chromium；nanowells：Wafergen iCELL8, SeqW ell8, CelSee）的发展也使得单细胞测序费用更加亲民。

然而，目前在大规模数据集分析中的主要挑战是将肿瘤细胞在基质和免疫细胞区分出来独立研究。

先前提出的inferCNV（测试过）和HoneyBadger（未测试）的方法，可以在具有足够大的基因组覆盖度的数据中估计基因组CNV，其更适用于第一代单细胞测序技术（在此理解为这两种方法更适用于长读长、测序深度足够的单细胞测序技术），而对于仅对mRNA的3 '或5 '端进行较低深度的高通量单细胞转录组技术检测效力不太适合[1]。

在2020年9月份发表的关于鼻咽癌单细胞的研究中，也可看出利用inferCNV工具分别对smart-seq2的和10x数据进行CNV检测，smart-seq2的信号要强于10x。

另外，先前方法不能准确地预测染色体断点的基因组位置，也不能根据肿瘤细胞和正常细胞的非整倍体拷贝数对细胞进行分类。在使用inferCNV时可能需要我们认为划定阈值，并且结合肿瘤样本和正常样本的比例，来确定肿瘤细胞和正常细胞。

在这里，给大家介绍一款软件—CopyKat，来区分肿瘤的拷贝数变异。

CopyKat 原 理 和 流 程

CopyKat的算法原理和流程大概理解一下：

• 将单细胞转录组的基因表达矩阵作为输入，依据基因组坐标对基因进行注释、排序；利用Freeman–Tukey transformation来稳定方差；使用polynomial dynamic linear modeling (DLM)矫正矩阵中异常值。

• 为了获得高置信度的二倍体细胞（正常细胞），先将细胞分成几个小的不同层次的亚群中，并利用 Gaussian mixture model (GMM)估计每个层次亚群的方差，具有最小估计方差的亚群被认为是高置信度的二倍体细胞。

• 为了检测染色体断点，用Poisson-gamma模型和Markov chain Monte Carlo (MCMC)迭代计算每个基因窗口的后验均值，并用Kolmogorov–Smirnov(KS)检验将没有显著性差异的window合并。另，为了加速计算，先将细胞分成不同群，检测共同的染色体断点并合并形成细胞群的基因组断点。

• 通过对单细胞拷贝数矩阵聚类，将细胞分成的不同亚群，并计算非整倍体细胞（肿瘤细胞）和二倍体基质细胞（正常细胞）之间的最大距离，对于距离不显著的需要使用GMM模型预测单个肿瘤细胞。

• 将单细胞拷贝数数据聚类以鉴定克隆亚群，可进一步分析不同肿瘤亚群的基因表达差异。

CopyKat 使 用

1、安装

library(devtools)
install_github("navinlabcode/copykat")

对于在集群上安装可能会存在网络限制以及权限限制，可以按照下面的方法尝试：

1、在Github网站下载copykat的zip文件；

2、放到集群上，解压缩，利用R CMD build folder 将其编译成二进制的包；

3、生成的tar.gz的R包就可以通过相应的R版本本地安装。

2、准备表达矩阵作为输入

 library(Seurat)
 raw <- Read10X(data.dir = data.path.to.cellranger.outs)
 raw <- CreateSeuratObject(counts = raw, project = "copycat.test", min.cells = 0, min.features = 0)
 exp.rawdata <- as.matrix(raw@assays$RNA@counts)
 write.table(exp.rawdata, file="exp.rawdata.txt", sep="\t", quote = FALSE, row.names = TRUE)

在这，我用的是RDS文件中提取的数据，省去重新创建Seurat对象的步骤。

immune.combined<-readRDS("immune.combined_umap.rds")
counts_matrix = as.matrix(immune.combined@assays$RNA@counts)
write.table(counts_matrix,file="counts_matrix_new.txt",col.names=T,row.names=T,sep="\t",quote=F)

3、运行copykat程序

export R_LIBS="/DL_and_ML/library/" && ymcR
library(copykat)
copykat.test <- copykat(rawmat=counts_matrix, id.type="S", cell.line="no", ngene.chr=5, win.size=25, KS.cut=0.2, sam.name="test", distance="euclidean", n.cores=10) ##运行速度1675个细胞六分钟左右

4、预览预测结果

第一部分结果是对细胞进行正常细胞和肿瘤细胞的判定结果，如下图。

pred.test <- data.frame(copykat.test$prediction)
head(pred.test)

第二部分结果是每个CNV区段在不同细胞的表达矩阵，前三列是segment在基因组上的位置信息。

CNA.test <- data.frame(copykat.test$CNAmat)
head(CNA.test[ , 1:5])
#Copykat可以生成用于衡量CNV水平的热图，行为细胞，列为基因组位置
my_palette <- colorRampPalette(rev(RColorBrewer::brewer.pal(n = 3, name = "RdBu")))(n = 999)
chr <- as.numeric(CNA.test$chrom) %% 2+1
rbPal1 <- colorRampPalette(c('grey60','grey30'))
CHR <- rbPal1(2)[as.numeric(chr)]
chr1 <- cbind(CHR,CHR)
rbPal5 <- colorRampPalette(RColorBrewer::brewer.pal(n = 8, name = "Dark2")[2:1])
com.preN <- pred.test$copykat.pred
pred <- rbPal5(2)[as.numeric(factor(com.preN))]

cells <- rbind(pred,pred)
col_breaks=c(seq(-1,-0.4,length=50),seq(-0.4,-0.2,length=150),seq(-0.2,0.2,length=600),seq(0.2,0.4,length=150),seq(0.4, 1,length=50))
pdf("CNV_heatmap.pdf",8,10)
heatmap.3(t(CNA.test[,4:ncol(CNA.test)]),dendrogram="r", distfun = function(x) parallelDist::parDist(x,threads =4, method = "euclidean"), hclustfun = function(x) hclust(x, method="ward.D2"), ColSideColors=chr1,RowSideColors=cells,Colv=NA, Rowv=TRUE,
 notecol="black",col=my_palette,breaks=col_breaks, key=TRUE,keysize=1, density.info="none", trace="none",
 cexRow=0.1,cexCol=0.1,cex.main=1,cex.lab=0.1, symm=F,symkey=F,symbreaks=T,cex=1, cex.main=4, margins=c(10,10))
legend("topright", paste("pred.",names(table(com.preN)),sep=""), pch=15,col=RColorBrewer::brewer.pal(n = 8, name = "Dark2")[2:1], cex=0.6, bty="n")
dev.off()

5、定义非整倍体肿瘤细胞的亚群

由上面的步骤可观察到正常细胞（diploid）和肿瘤细胞（aneuploid），接下来可以提取判定为肿瘤细胞的所有细胞，通过聚类将肿瘤细胞再分成不同的亚群，由此来研究不同肿瘤细胞亚群的差异。

tumor.cells <- pred.test$cell.names[which(pred.test$copykat.pred=="aneuploid")]
tumor.mat <- CNA.test[, which(colnames(CNA.test) %in% tumor.cells)]
hcc <- hclust(parallelDist::parDist(t(tumor.mat),threads =4, method = "euclidean"), method = "ward.D2")
hc.umap <- cutree(hcc,2)
rbPal6 <- colorRampPalette(RColorBrewer::brewer.pal(n = 8, name = "Dark2")[3:4])
subpop <- rbPal6(2)[as.numeric(factor(hc.umap))]
cells <- rbind(subpop,subpop)
pdf("Subpopulation_CNV_heatmap.pdf",8,10)
heatmap.3(t(tumor.mat),dendrogram="r", distfun = function(x) parallelDist::parDist(x,threads =4, method = "euclidean"), hclustfun = function(x) hclust(x, method="ward.D2"), ColSideColors=chr1,RowSideColors=cells,Colv=NA, Rowv=TRUE,notecol="black",col=my_palette,breaks=col_breaks, key=TRUE, keysize=1, density.info="none", trace="none", cexRow=0.1,cexCol=0.1,cex.main=1,cex.lab=0.1, symm=F,symkey=F,symbreaks=T,cex=1, cex.main=4, margins=c(10,10))
legend("topright", c("c1","c2"), pch=15,col=RColorBrewer::brewer.pal(n = 8, name = "Dark2")[3:4], cex=0.9, bty='n')
dev.off()

这里用的不是项目分析中的实际数据，就不做批次效应的矫正了。

rds <- immune.combined
rds <- NormalizeData(rds,verbose=F)
rds <- ScaleData(rds,verbose=F)
rds <- FindVariableFeatures(rds,verbose=F)
rds <- RunPCA(rds, npcs = 30, verbose = F)
rds_umap <- RunUMAP(rds, reduction = "pca",dims = 1:30)
rds_umap <- FindNeighbors(object = rds_umap , reduction = "pca", dims = 1:30)
rds_umap <- FindClusters(rds_umap , resolution = 0.5)

在运行正常/肿瘤细胞预测时可能是做了过滤，测试项目的细胞共1675个，而预测出来只有1672个。经查看，缺少了下面三个细胞。

rownames([e-mail])[!rownames([e-mail]) %in% pred.test$cell.names]
\#[1] "GTGTGATAGGCCATAG-1" "GCTTGGGAGGCTAACG-1" "GTACAGTAGTCCTGCG-1"

防止绘制聚类图时出错，我们先对pred.test文件将缺失的三个细胞补上。

omit_cell <- data.frame(cell.names=c("GTGTGATAGGCCATAG-1","GCTTGGGAGGCTAACG-1","GTACAGTAGTCCTGCG-1"),copykat.pred=rep("unknown",3))
rownames(omit_cell)<-c("GTGTGATAGGCCATAG-1","GCTTGGGAGGCTAACG-1","GTACAGTAGTCCTGCG-1")
pred.test1 <-rbind(pred.test,omit_cell)
pred.test1$cell.names<-gsub('-','.',pred.test1$cell.names)
rownames(pred.test1)<-gsub('-','.',rownames(pred.test1))


[e-mail]$cell.names0<-rownames([e-mail])
[e-mail]$cell.names<-gsub('-','.',rownames([e-mail]))
[e-mail]<-merge([e-mail],pred.test1,by="cell.names")
rownames([e-mail])<[e-mail]$cell.names0
[e-mail]$tumor_subpopulation <- "NA"
normal.cells <- gsub('-','.',pred.test$cell.names)[which(pred.test$copykat.pred=="diploid")]
[e-mail]$tumor_subpopulation[[e-mail]$cell.names %in% normal.cells] <- "normal"
[e-mail]$tumor_subpopulation[[e-mail]$cell.names %in% names(hc.umap[hc.umap==1])] <- "tumor_subpopulation1"
[e-mail]$tumor_subpopulation[[e-mail]$cell.names %in% names(hc.umap[hc.umap==2])] <- "tumor_subpopulation2"

color_all <-c("#4DBBD5FF","#E64B35FF","#00A087FF","#3C5488FF","#F39B7FFF","#8491B4FF","#91D1C2FF","#DC0000FF","#7E6148FF","#B09C85FF","#3B4992FF","#EE0000FF","#008B45FF")
length_group <-length(unique([e-mail]$seurat_clusters))
length_pred <- length(unique([e-mail]$copykat.pred))
length_subpopulation<- length(unique([e-mail]$tumor_subpopulation))


p0<-theme(panel.grid=element_blank(), legend.background = element_rect(colour = NA),
 legend.title = element_blank(),legend.text = element_text(face="plain", color="black",size=20),
 axis.text.x = element_text(color="black",size=20),
 axis.text.y = element_text(color="black",size=20),
 axis.title.x = element_text(face="plain", color="black",size=20),
 axis.title.y = element_text(face="plain", color="black",size=20))
pdf("umap_cluster_samples.pdf",20,12)

p1 <- DimPlot(rds_umap, reduction = "umap", cols=color_all[1:length_group],label = TRUE,label.size = 8,pt.size = 1)+p0
p2 <- DimPlot(rds_umap, reduction = "umap", cols=color_all[1:length_pred],group.by = "copykat.pred",pt.size = 1)+p0
p3 <- DimPlot(rds_umap, reduction = "umap", cols=color_all[1:length_subpopulation],group.by = "tumor_subpopulation",pt.size = 1)+p0
p4<-plot_grid(p1,p2,p3,ncol=3)
print(p4)
dev.off()

总  结

亲测过inferCNV和CopyKAT，在计算上CopyKAT是基于基因组的坐标位置来检测CNV水平，而是基因的表达量，对于计算区间的CNV水平相较于单个基因来说，噪音基因的影响更小。

在计算速度上，CopyKAT可以利用并行提高运算速度，虽未拿同一数据比较，但从等待的焦虑程度来说，CopyKAT等待的焦虑程度更低一些。

总之，值得一试。

参  考

1. https://github.com/navinlabcode/copykat

2. Gao, R., Bai, S., Henderson, Y. C., Lin, Y., Schalck, A., Yan, Y., Kumar, T., Hu, M., Sei, E., Davis, A., Wang, F., Shaitelman, S. F., Wang, J. R., Chen, K., Moulder, S., Lai, S. Y. & Navin, N. E. (2021). Delineating copy number and clonal substructure in human tumors from single-cell transcriptomes. Nat Biotechnol. doi:10.1038/s41587-020-00795-2.

作者：Bio_gevin

审稿：童蒙

编辑：angelica

往 期 回 顾